Set up

Loading packages

knitr::opts_chunk$set(echo = TRUE)
#install.packages("pacman")
pacman::p_load(here, dplyr, ggplot2, ggpubr, viridis, ggh4x, stringr,
               forcats, tidyr, lme4, emmeans, plotly)
# Ryssa's go-to theme
theme_set(theme_light()+
            theme(
              plot.title = element_text(size=rel(1.2),face="bold"),
              axis.title = element_text(size=rel(1),face="bold"),
              axis.text = element_text(size=rel(1),colour = 'black'),
              strip.text = element_text(size=rel(1),colour = 'black', 
                                        face = "bold"),
              legend.text = element_text(size=rel(1)),
              legend.title = element_text(size=rel(1),face="bold"),
              panel.grid = element_blank()))

Data sets

Movement synchrony

Pairs of people played the mirror game. We measured how well people were able to mirror each other’s spontaneous movements to understand how movement synchrony influences response inhibition and the related patterns of brain activity.

Moffat et al. (2024). Inhibiting responses under the watch of a recently synchronized peer increases self-monitoring: evidence from functional near-infrared spectroscopy.

For this example, we keep the columns:

  • ID participant identifier
  • leader whether the participant (P) or helper (C) was leading the game
  • meanSim the mean level of movement synchrony during the game
  • group whether the participant belonged to the synchronised or control group
poseSimilarity <- read.csv(here::here("data/pose_similarity.csv")) %>%
  filter(ID <500) %>%
  select(-c(cutoff, nFrameskip, nFrames, sdSim)) %>%
  mutate(group = factor(case_when(ID <200 ~ "Synchronised", TRUE ~ "Control")),
         leader = factor(case_when(leader == "C" ~ "Helper",
                                   leader == "P" ~ "Participant")))

head(poseSimilarity)
##    ID      leader   meanSim        group
## 1 101      Helper 0.7625684 Synchronised
## 2 101 Participant 0.8129776 Synchronised
## 3 102      Helper 0.7452378 Synchronised
## 4 102 Participant 0.8203161 Synchronised
## 5 103      Helper 0.8224377 Synchronised
## 6 103 Participant 0.7534514 Synchronised

Synchrony + Traits

While people played the mirror game, we tracked the level of synchrony and the complexity of the mirror game movements. The participants also filled in questionnaires about themselves. We were particularly interested in embodied movement knowledge. One such measure is “body competence”, which indexes how strongly people believe that their body can accomplish physical activities.

Moffat et al. (2024). Dyadic body competence predicts movement synchrony during the mirror game.

For this example, we keep the columns:

  • ID participant identifier
  • M_similarity mean movement synchrony
  • M_entropy mean movement complexity computed using entropy
  • mean_BCQ mean body competence score per dyad
  • Role whether the participant was leading the game (leader) or following the helper (following)

We also create a new column:

  • BCQ_half to categorise dyads as being above or below the median of mean_BCQ

Prepare data

poses_traits <- read.csv(here("data/mirrorring_embodiment.csv")) %>%
    select(ID, M_similarity, M_entropy, BCQ, Conf_BCQ, Role, video_code) 

# Calculate mean and z-scores for helper and participant BCQ
poses_traits1 <- poses_traits %>% 
  # Calculate the mean separately for each row (pair)
  rowwise() %>% 
  mutate(mean_BCQ = mean(c(BCQ, Conf_BCQ)), # BCQ = body competence
         BCQ_half = factor(case_when(mean_BCQ > 7.5 ~ "High",
                              TRUE ~ "Low"))) %>%
  select(ID, M_similarity, M_entropy, mean_BCQ, BCQ_half, Role) 

head(poses_traits1)
## # A tibble: 6 × 6
## # Rowwise: 
##      ID M_similarity M_entropy mean_BCQ BCQ_half Role     
##   <int>        <dbl>     <dbl>    <dbl> <fct>    <chr>    
## 1   114        0.888    0.0555      9.5 High     Following
## 2   114        0.907    0.149       9.5 High     Following
## 3   114        0.629    0.107       9.5 High     Following
## 4   114        0.841    0.0384      9.5 High     Following
## 5   114        0.735    0.0637      9.5 High     Following
## 6   114        0.557    0.0874      9.5 High     Following

1. Model to figure

Let’s compare the level of synchrony observed during the mirror game (group = Synchronised) vs. the control game (group = Control), and find out if it matters who was leading (Leader or Helper in leader).

head(poseSimilarity)
##    ID      leader   meanSim        group
## 1 101      Helper 0.7625684 Synchronised
## 2 101 Participant 0.8129776 Synchronised
## 3 102      Helper 0.7452378 Synchronised
## 4 102 Participant 0.8203161 Synchronised
## 5 103      Helper 0.8224377 Synchronised
## 6 103 Participant 0.7534514 Synchronised

Model with lme4

model <- lmer(meanSim ~ group*leader + (1|ID), data = poseSimilarity)
summary(model)
## Linear mixed model fit by REML ['lmerMod']
## Formula: meanSim ~ group * leader + (1 | ID)
##    Data: poseSimilarity
## 
## REML criterion at convergence: -219
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.46578 -0.39503 -0.02317  0.46849  2.29187 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.002627 0.05125 
##  Residual             0.005684 0.07539 
## Number of obs: 120, groups:  ID, 60
## 
## Fixed effects:
##                                     Estimate Std. Error t value
## (Intercept)                          0.46831    0.01664  28.137
## groupSynchronised                    0.31627    0.02354  13.436
## leaderParticipant                   -0.07133    0.01947  -3.664
## groupSynchronised:leaderParticipant  0.11021    0.02753   4.003
## 
## Correlation of Fixed Effects:
##             (Intr) grpSyn ldrPrt
## grpSynchrns -0.707              
## ledrPrtcpnt -0.585  0.413       
## grpSynchr:P  0.413 -0.585 -0.707

Estimates with emmeans

model %>% emmeans(specs = ~ group + leader)
##  group        leader      emmean     SE  df lower.CL upper.CL
##  Control      Helper       0.468 0.0166 105    0.435    0.501
##  Synchronised Helper       0.785 0.0166 105    0.752    0.818
##  Control      Participant  0.397 0.0166 105    0.364    0.430
##  Synchronised Participant  0.823 0.0166 105    0.790    0.856
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95
model_estimates <- model %>% emmeans(specs = ~ group + leader) %>%
  as.data.frame()

head(model_estimates)
##  group        leader         emmean         SE     df  lower.CL  upper.CL
##  Control      Helper      0.4683098 0.01664421 105.46 0.4353091 0.5013105
##  Synchronised Helper      0.7845749 0.01664421 105.46 0.7515742 0.8175756
##  Control      Participant 0.3969787 0.01664421 105.46 0.3639780 0.4299794
##  Synchronised Participant 0.8234531 0.01664421 105.46 0.7904524 0.8564538
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

Plot estimates

Geom_point ()

model_estimates %>%
  ggplot(aes(x = group, y = emmean, color = leader))+
  geom_point(size = 4, position = position_dodge(0.7))+
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                position = position_dodge(0.7), width = 0.1)+
  scale_color_manual(values = c("#89398a", "#398a89"), name = "Leader")+
  ylim(0,1) +
  labs(x = "Group", y = "Parameter estimate and 95% CI")

### Geom_bar+point()

estimates <- model_estimates %>%
  ggplot(aes(x = group, y = emmean, fill = leader))+
  geom_bar(stat = "identity", alpha = 0.9, position = position_dodge(.9))+
  geom_point(size = 4, position = position_dodge(.9))+
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                position = position_dodge(.9), width = 0.1)+
  scale_fill_manual(values = c("#89398a", "#398a89"), name = "Leader")+
  ylim(0,1) +
  labs(x = "Group", y = "Parameter estimate and 95% CI")+
  theme(legend.position = c(0.25, 0.85))
## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
estimates

Contrasts with emmeans

model %>% emmeans(specs = ~ group) %>% pairs()
## NOTE: Results may be misleading due to involvement in interactions
##  contrast               estimate     SE df t.ratio p.value
##  Control - Synchronised   -0.371 0.0191 58 -19.449  <.0001
## 
## Results are averaged over the levels of: leader 
## Degrees-of-freedom method: kenward-roger
group_contrast <- model %>% emmeans(specs = ~ group) %>% 
  pairs() %>%
  as.data.frame()
## NOTE: Results may be misleading due to involvement in interactions
model %>% emmeans(specs = ~ leader) %>% pairs()
## NOTE: Results may be misleading due to involvement in interactions
##  contrast             estimate     SE df t.ratio p.value
##  Helper - Participant   0.0162 0.0138 58   1.179  0.2433
## 
## Results are averaged over the levels of: group 
## Degrees-of-freedom method: kenward-roger
leader_contrast <- model %>% emmeans(specs = ~ leader) %>%
  pairs() %>%
  as.data.frame()
## NOTE: Results may be misleading due to involvement in interactions
contrast_df <- rbind(group_contrast, leader_contrast) %>%
  mutate(contrast = case_when(contrast == "Helper - Participant" ~ "Helper\n-Participant", 
                              contrast == "Control - Synchronised" ~ "Control\n-Synchronised"))

head(contrast_df)
##                 contrast    estimate        SE df    t.ratio      p.value
## 1 Control\n-Synchronised -0.37136976 0.0190943 58 -19.449245 4.311649e-27
## 2   Helper\n-Participant  0.01622648 0.0137647 58   1.178847 2.432725e-01

Plot contrasts

contrasts <- contrast_df %>%
  ggplot(aes(x = contrast, y = estimate, fill = contrast))+
  geom_hline(yintercept = 0, color = "grey", linetype = "dashed")+
  geom_errorbar(aes(ymin = estimate-(SE*1.96), #*1.96 for condidence intervals
                    ymax = estimate+(SE*1.96), width = 0.1))+  
  geom_point(size = 4, shape =21)+
  scale_fill_manual(values = c("magenta", "#80ED99"))+
  labs(x = "Contrast estimate", y = "Contrast")+
  theme(legend.position = "none")+
  annotate(x = 1, y = -0.3, label = "***", geom = "text")

contrasts

Combine plots

In the rendered markdown, the labels overlap. This is not the case in the saved image.

ggarrange(estimates, contrasts)

ggsave("figures/modelplots.jpg", height = 4, width = 9, dpi = 400)

2. 3D plots

3D plots can help visualise data with multiple continuous dimensions. The data set with movement synchrony, movement complexity and dyad-level measures of body competence is perfect for this example.

head(poses_traits1)
## # A tibble: 6 × 6
## # Rowwise: 
##      ID M_similarity M_entropy mean_BCQ BCQ_half Role     
##   <int>        <dbl>     <dbl>    <dbl> <fct>    <chr>    
## 1   114        0.888    0.0555      9.5 High     Following
## 2   114        0.907    0.149       9.5 High     Following
## 3   114        0.629    0.107       9.5 High     Following
## 4   114        0.841    0.0384      9.5 High     Following
## 5   114        0.735    0.0637      9.5 High     Following
## 6   114        0.557    0.0874      9.5 High     Following

Continuous color scale

The continuous color scale for one axis can help give the plot depth, and thereby help interpretation when rotating it. Click on the plot and drag it around to see different orientations of the axes.

poses_traits1 %>%
  plot_ly(x = ~M_similarity, y = ~mean_BCQ, z = ~M_entropy,
          marker = list(color = ~M_similarity,
                        colorscale = c('#FFE1A1', '#683531'),
                        showscale = TRUE)) %>%
  add_markers() %>%
  layout(scene = list(xaxis = list(title = 'Synchrony'),
                                   yaxis = list(title = 'Dyadic Body Competence'),
                                   zaxis = list(title = 'Complexity')),
                      annotations = list(
                        x = 1.13,
                        y = 1.05,
                        text = 'Synchrony',
                        xref = 'paper',
                        yref = 'paper',
                        showarrow = FALSE
                        ))

Discrete colors

Alternatively, by creating a new column with which to divide the data into categories, we can help out interpretation too.Click on the plot and drag it around to see different orientations of the axes.

poses_traits1 %>%
  plot_ly(x = ~M_similarity, y = ~mean_BCQ, z = ~M_entropy,
          color = ~BCQ_half, colors = c("#89398a", "#398a89"), opacity = 0.5) %>%
  add_markers() %>%
  layout(scene = list(xaxis = list(title = 'Synchrony'),
                     yaxis = list(title = 'Dyadic Body Compentence'),
                     zaxis = list(title = 'Complexity')))